{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script type=\"text/x-mathjax-config\">\n",
    "  MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "  });\n",
    "  MathJax.Hub.Config({\n",
    "    TeX: { extensions: [\"AMSmath.js\", \"AMSsymbols.js\", \"autobold.js\", \"autoload-all.js\"] }\n",
    "  });\n",
    "  MathJax.Hub.Config({\n",
    "    tex2jax: {\n",
    "      inlineMath: [['$','$'], ['\\\\(','\\\\)']],\n",
    "      processEscapes: true\n",
    "    }\n",
    "  });\n",
    "</script>\n",
    "<script type=\"text/javascript\" src=\"http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML\">\n",
    "</script>\n",
    "\n",
    "# User guide\n",
    "\n",
    "---\n",
    "\n",
    "`TaylorSeries.jl` can be thought of as a polynomial algebraic manipulator in one or more\n",
    "variables; these two cases are treated separately for now.  Three new types are defined,\n",
    "`Taylor1`, `HomogeneousPolynomial` and `TaylorN`, which correspond to\n",
    "expansions in one independent variable, homogeneous polynomials of various variables, and the polynomial\n",
    "series in many independent variables, respectively. These types are subtypes\n",
    "of `Number` and are defined parametrically.\n",
    "\n",
    "The package is loaded as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: New definition \n",
      "    /(TaylorSeries.Taylor1{#T<:Real}, #T<:Real) at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:261\n",
      "is ambiguous with: \n",
      "    /(TaylorSeries.Taylor1{Base.Rational{#T<:Integer}}, #S<:Union{Base.Complex, Real}) at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:254.\n",
      "To fix, define \n",
      "    /(TaylorSeries.Taylor1{_<:Base.Rational{#T<:Integer}}, _<:Base.Rational{#T<:Integer})\n",
      "before the new definition.\n"
     ]
    }
   ],
   "source": [
    "using TaylorSeries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One variable\n",
    "\n",
    "Taylor expansions in one variable are represented by the `Taylor1` type, which\n",
    "consists of a vector of coefficients (field `coeffs`) and the maximum\n",
    "order considered for the expansion (field `order`). The\n",
    "coefficients are arranged in ascending order with respect to the power of the\n",
    "independent variable, so that\n",
    "`coeffs[1]` is the constant term, `coeffs[2]` gives the first order term,\n",
    "etc. This is a dense representation of the polynomial.\n",
    "The order of the polynomial can be\n",
    "omitted in the constructor, which is then fixed from the length of the\n",
    "vector of coefficients; otherwise, the maximum\n",
    "of the length of the vector of coefficients and the given integer is taken."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1 + 2 t + 3 t² + 𝒪(t³)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor1([1, 2, 3]) # Polynomial of order 2 with coefficients 1, 2, 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " ( 1.0 im ) t + 𝒪(t²)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor1([0.0, 1im]) # Also works with complex numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "affine (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "affine(a) = a + Taylor1(typeof(a),5)  ## a + t of order 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 t + 𝒪(t⁶)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = affine(0.0) # Independent variable `t`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the information about the maximum order considered is displayed\n",
    "using a big-O notation.\n",
    "\n",
    "The definition of `affine(a)` uses the method `Taylor1{T<:Number}(::Type{T},order::Int64)`, which is a\n",
    "shortcut to define the independent variable of a Taylor expansion,\n",
    "with a given type and given order. As we show below, this is one of the\n",
    "easiest ways to work with the package.\n",
    "\n",
    "The usual arithmetic operators (`+`, `-`, `*`, `/`, `^`, `==`) have been\n",
    "extended to work with the `Taylor1` type, including promotions that involve\n",
    "`Number`s. The operations return a valid Taylor expansion with the same\n",
    "maximum order; compare the last example below, where this is not possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 2.5 t + 3.0 t² + 𝒪(t⁶)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t*(3t+2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 t + 1.0 t² + 1.0 t³ + 1.0 t⁴ + 1.0 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/(1-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " - 2.0 t + 1.0 t² + 𝒪(t⁶)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t*(t^2-4)/(t+2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " ( 1.0 im ) t + 𝒪(t⁶)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tI = im*t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 0.0 + 𝒪(t⁶)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t^6  # order is 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 - 3.2 t + 3.5200000000000005 t² - 1.4080000000000004 t³ + 0.07040000000000009 t⁴ + 0.011264000000000012 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1-t)^3.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 t² - 0.5 t³ + 0.8333333333333333 t⁴ - 0.75 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1+t)^t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If no valid Taylor expansion can be computed, an error is thrown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ArgumentError: Division does not define a Taylor1 polynomial\nor its first non-zero coefficient is Inf/NaN.\nOrder k=0 => coeff[1]=Inf.",
     "output_type": "error",
     "traceback": [
      "ArgumentError: Division does not define a Taylor1 polynomial\nor its first non-zero coefficient is Inf/NaN.\nOrder k=0 => coeff[1]=Inf.",
      "",
      " in divfactorization at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:308",
      " in / at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:284",
      " in / at promotion.jl:170"
     ]
    }
   ],
   "source": [
    "1/t "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ArgumentError: The 0th order Taylor1 coefficient must be non-zero\nto raise the Taylor1 polynomial to a non-integer exponent.",
     "output_type": "error",
     "traceback": [
      "ArgumentError: The 0th order Taylor1 coefficient must be non-zero\nto raise the Taylor1 polynomial to a non-integer exponent.",
      "",
      " in ^ at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:450"
     ]
    }
   ],
   "source": [
    "t^3.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several elementary functions have been implemented; these compute their\n",
    "coefficients recursively. So far, these functions are `exp`, `log`, `sqrt`, `sin`, `cos`\n",
    "and `tan`;\n",
    "more will be added in the future. Note that this way of obtaining the\n",
    "Taylor coefficients is not the *laziest* way, in particular for many independent\n",
    "variables. Yet, it is quite efficient, especially for the integration of\n",
    "ordinary differential equations, which is among the applications we have in mind."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " - 1.0 t - 0.5 t² - 0.3333333333333333 t³ - 0.25 t⁴ - 0.2 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log(1-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ArgumentError: First non-vanishing Taylor1 coefficient must correspond\nto an **even power** in order to expand `sqrt` around 0.",
     "output_type": "error",
     "traceback": [
      "ArgumentError: First non-vanishing Taylor1 coefficient must correspond\nto an **even power** in order to expand `sqrt` around 0.",
      "",
      " in sqrt at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:563"
     ]
    }
   ],
   "source": [
    "sqrt(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 0.5 t - 0.125 t² + 0.0625 t³ - 0.0390625 t⁴ + 0.02734375 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqrt(1 + t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " - 1.0 t + 0.16666666666666666 t³ - 0.008333333333333333 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imag(exp(tI)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "real(exp(Taylor1([0.0,1im],17))) - cos(Taylor1([0.0,1.0],17)) == 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1//1 + 1//1 t + 1//2 t² + 1//6 t³ + 1//24 t⁴ + 1//120 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Taylor1{Rational{Int64}}, exp(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Differentiating and integrating is straightforward for polynomial expansions in\n",
    "one variable. The last coefficient of a derivative is set to zero to keep the\n",
    "same order as the original polynomial; for the integral, an\n",
    "integration constant may be set to a different value (the default is zero). The\n",
    "order of the resulting polynomial is not changed.\n",
    "The $n$-th ($n \\ge 0$) derivative at `t=0` is obtained using `derivative(n,a)`  where `a` is a Taylor series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "derivative(exp(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(exp(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate( exp(t), 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate( derivative( exp(-t)), 1.0 ) == exp(-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate(derivative(exp(affine(1.0))),0.) == exp(1.0)  # deriv of `exp(1+t)` at t=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "derivative(5, exp(affine(1.0))) == exp(1.0) # Fifth derivative of `exp(1+t) at t=0`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate a Taylor series at a point, Horner's rule is used via the function\n",
    "`evaluate(a::Taylor, dt::Number)`. Here, $dt$ is the increment from\n",
    "the point $t_0$ where the Taylor expansion is calculated, i.e., the series\n",
    "is evaluated at $t = t_0 + dt$. Omitting $dt$ corresponds to $dt = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate(exp(affine(1.0))) - e # exp(t) around t0=1 (order 5), evaluated there (dt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0016151617923783057"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate(exp(t), 1) - e # exp(t) around t0=0 (order 5), evaluated at t=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate( exp( Taylor1(Float64,17) ), 1) - e # exp(t) around t0=0 (order 17), evaluated at t=1\n",
    "0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.000000000000000000000000000000000000000000000000000000000000000000000000000000 t + 𝒪(t⁵¹)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tBig = Taylor1([zero(BigFloat),one(BigFloat)],50) # With BigFloats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.718281828459045235360287471352662497757247093699959574966967627723419298053556"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eBig = evaluate( exp(tBig), one(BigFloat) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.573322999985292556154129119543257102601105719980995128942636339920549561322098e-67"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e - eBig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Many variables\n",
    "\n",
    "A polynomial in $N>1$ variables can be represented in (at least) two ways:\n",
    "As a vector whose coefficients are homogeneous polynomials of fixed degree, or\n",
    "as a vector whose coefficients are polynomials in $N-1$ variables. We have opted\n",
    "to implement the first option, which seems to show better performance. An elegant\n",
    "(lazy) implementation of the second representation was discussed on the\n",
    "[julia-users](https://groups.google.com/forum/#!msg/julia-users/AkK_UdST3Ig/sNrtyRJHK0AJ) list.\n",
    "\n",
    "`TaylorN` is thus constructed as a vector of parameterized homogeneous polynomials\n",
    "defined by the type `HomogeneousPolynomial`, which in turn is a vector of\n",
    "coefficients of given order (degree). This implementation imposes that the user\n",
    "has to specify the (maximum) order and the number of independent\n",
    "variables, which is done using the `set_variables(names)` function.\n",
    "`names` is a string consisting of the desired *output* names of the variables,\n",
    "separated by spaces. A vector of the resulting Taylor variables is returned:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{TaylorSeries.TaylorN{Float64},1}:\n",
       "  1.0 x + 𝒪(‖x‖⁷)\n",
       "  1.0 y + 𝒪(‖x‖⁷)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = set_variables(\"x y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 x + 𝒪(‖x‖⁷)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TaylorSeries.TaylorN{Float64}"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7-element Array{TaylorSeries.HomogeneousPolynomial{Float64},1}:\n",
       "    0.0\n",
       "  1.0 x\n",
       "    0.0\n",
       "    0.0\n",
       "    0.0\n",
       "    0.0\n",
       "    0.0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown, the resulting objects are of `TaylorN{Float64}` type. There is an optional `order` keyword argument for `set_variables`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{TaylorSeries.TaylorN{Float64},1}:\n",
       "  1.0 x + 𝒪(‖x‖¹¹)\n",
       "  1.0 y + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set_variables(\"x y\", order=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numbered variables are also available by specifying a single\n",
    "variable name and the optional keyword argument `numvars`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{TaylorSeries.TaylorN{Float64},1}:\n",
       "  1.0 α₁ + 𝒪(‖x‖⁷)\n",
       "  1.0 α₂ + 𝒪(‖x‖⁷)\n",
       "  1.0 α₃ + 𝒪(‖x‖⁷)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set_variables(\"α\", numvars=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``show_params_TaylorN()` displays the current values of the parameters, in an info block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[34mINFO: Parameters for `TaylorN` and `HomogeneousPolynomial`:\n",
      "Maximum order       = 6\n",
      "Number of variables = 3\n",
      "Variable names      = Union{ASCIIString,UTF8String}[\"α₁\",\"α₂\",\"α₃\"]\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "show_params_TaylorN()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Internally, changing these parameters defines dictionaries that\n",
    "translate the position of the coefficients of a `HomogeneousPolynomial`\n",
    "into the corresponding\n",
    "multi-variable monomials. Fixing these values from the start is imperative.\n",
    "\n",
    "The easiest way to construct a `TaylorN` object is by defining symbols for\n",
    "the independent variables, as above. Again, the Taylor expansions are implemented\n",
    "around 0 for all variables; if the expansion\n",
    "is needed around a different value, the trick is a simple translation of\n",
    "the corresponding\n",
    "independent variable $x \\to x+a$.\n",
    "\n",
    "\n",
    "Other ways of constructing `TaylorN` polynomials involve using `HomogeneousPolynomial`\n",
    "objects directly, which is uncomfortable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set_variables(\"x\", numvars=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1 x₁ - 1 x₂"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HomogeneousPolynomial([1,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1 x₁ + 1 x₁² + 2 x₁ x₂ + 3 x₂² + 𝒪(‖x‖⁵)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TaylorN( [HomogeneousPolynomial([1,0]), HomogeneousPolynomial([1,2,3])], 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, the usual arithmetic operators (`+`, `-`, `*`, `/`, `^`, `==`)\n",
    "have been extended to work with `TaylorN` objects, including the appropriate\n",
    "promotions to deal with numbers. (Some of the arithmetic operations have\n",
    "also been extended for\n",
    "`HomogeneousPolynomial`, whenever the result is a `HomogeneousPolynomial`;\n",
    "division, for instance, is not extended.) Also, the elementary functions have been\n",
    "implemented, again by computing their coefficients recursively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, y = set_variables(\"x y\", order=10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 1.0 x + 1.0 y + 0.5 x² + 1.0 x y + 0.5 y² + 0.16666666666666666 x³ + 0.5 x² y + 0.5 x y² + 0.16666666666666666 y³ + 0.041666666666666664 x⁴ + 0.16666666666666666 x³ y + 0.25 x² y² + 0.16666666666666666 x y³ + 0.041666666666666664 y⁴ + 0.008333333333333333 x⁵ + 0.041666666666666664 x⁴ y + 0.08333333333333333 x³ y² + 0.08333333333333333 x² y³ + 0.041666666666666664 x y⁴ + 0.008333333333333333 y⁵ + 0.0013888888888888887 x⁶ + 0.008333333333333331 x⁵ y + 0.020833333333333332 x⁴ y² + 0.027777777777777776 x³ y³ + 0.020833333333333332 x² y⁴ + 0.008333333333333331 x y⁵ + 0.0013888888888888887 y⁶ + 0.00019841269841269839 x⁷ + 0.0013888888888888885 x⁶ y + 0.004166666666666666 x⁵ y² + 0.006944444444444443 x⁴ y³ + 0.006944444444444443 x³ y⁴ + 0.004166666666666666 x² y⁵ + 0.0013888888888888885 x y⁶ + 0.00019841269841269839 y⁷ + 2.4801587301587298e-5 x⁸ + 0.00019841269841269836 x⁷ y + 0.0006944444444444443 x⁶ y² + 0.0013888888888888887 x⁵ y³ + 0.0017361111111111108 x⁴ y⁴ + 0.0013888888888888887 x³ y⁵ + 0.0006944444444444443 x² y⁶ + 0.00019841269841269836 x y⁷ + 2.4801587301587298e-5 y⁸ + 2.7557319223985884e-6 x⁹ + 2.4801587301587295e-5 x⁸ y + 9.920634920634918e-5 x⁷ y² + 0.0002314814814814814 x⁶ y³ + 0.0003472222222222221 x⁵ y⁴ + 0.0003472222222222221 x⁴ y⁵ + 0.0002314814814814814 x³ y⁶ + 9.920634920634918e-5 x² y⁷ + 2.4801587301587295e-5 x y⁸ + 2.7557319223985884e-6 y⁹ + 2.7557319223985883e-7 x¹⁰ + 2.7557319223985884e-6 x⁹ y + 1.2400793650793647e-5 x⁸ y² + 3.306878306878306e-5 x⁷ y³ + 5.787037037037036e-5 x⁶ y⁴ + 6.944444444444443e-5 x⁵ y⁵ + 5.787037037037036e-5 x⁴ y⁶ + 3.306878306878306e-5 x³ y⁷ + 1.2400793650793647e-5 x² y⁸ + 2.7557319223985884e-6 x y⁹ + 2.7557319223985883e-7 y¹⁰ + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exy = exp(x+y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `getcoeff(a,v)`\n",
    "gives the coefficient of `a::TaylorN` that corresponds to the monomial\n",
    "specified by the vector of powers `v`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TaylorSeries.TaylorN{Float64}"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(exy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1//720"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rationalize(getcoeff(exy, [3,5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Partial differentiation is also implemented for `TaylorN` objects,\n",
    "through `derivative`; integration is yet to be implemented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f (generic function with 1 method)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x,y) = x^3 + 2x^2 * y - 7x + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "g (generic function with 1 method)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g(x,y) = y - x^4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " - 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "derivative( f(x,y), 1 )   # partial derivative with respect to 1st variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "derivative( g(x,y), 2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "AssertionError: 1 <= r <= get_numvars()",
     "output_type": "error",
     "traceback": [
      "AssertionError: 1 <= r <= get_numvars()",
      "",
      " in derivative at /Users/benet/.julia/v0.4/TaylorSeries/src/TaylorN.jl:948",
      " in derivative at /Users/benet/.julia/v0.4/TaylorSeries/src/TaylorN.jl:982"
     ]
    }
   ],
   "source": [
    "derivative( g(x,y), 3 )   # error, since we are dealing with 2 variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`evaluate` can also be used for `TaylorN` objects, using it on vectors of\n",
    "numbers (`Real` or `Complex`); the length of the vector must coincide with the number\n",
    "of independent variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate(exy, [.1,.02]) == e^0.12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions to compute the gradient, Jacobian and\n",
    "Hessian have also been implemented. Using the\n",
    "functions $f(x,y) = x^3 + 2x^2 y - 7 x + 2$ and $g(x,y) = y-x^4$ defined above,\n",
    "we may use `gradient` or `∇` (`\\nabla+TAB`); the results are of\n",
    "type `Array{TaylorN{T},1}`. To compute the Jacobian or Hessian of a vector field\n",
    "evaluated at a point, we use `jacobian` and `hessian`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 2.0 - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f1 = f(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 y - 1.0 x⁴ + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g1 = g(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{TaylorSeries.TaylorN{Float64},1}:\n",
       "  - 4.0 x³ + 𝒪(‖x‖¹¹)\n",
       "       1.0 + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient( g1 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{TaylorSeries.TaylorN{Float64},1}:\n",
       "  - 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)\n",
       "                    2.0 x² + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "∇(f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 2.0 - 7.0 x - 5.0 y + 14.0 x y + 1.0 x³ + 2.0 x² y + 5.0 x⁴ - 2.0 x³ y - 4.0 x² y² - 14.0 x⁵ + 2.0 x⁷ + 4.0 x⁶ y + 𝒪(‖x‖¹¹)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fg = f1-g1-2*f1*g1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2x2 Array{Float64,2}:\n",
       "  13.0  8.0\n",
       " -32.0  1.0"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobian([f1,g1], [2,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2x2 Array{Float64,2}:\n",
       " -26.0  20.0\n",
       "  20.0  -8.0"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hessian(fg, [1.0,1.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "### 1. Four-square identity\n",
    "\n",
    "Euler proved the [following four-square identity](https://en.wikipedia.org/wiki/Euler%27s_four-square_identity):\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "(a_1^2+a_2^2+a_3^2+a_4^2) \\cdot (b_1^2+b_2^2+b_3^2+b_4^2) = & &\n",
    "   (a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4)^2  \\\\\n",
    "  &+& (a_1 b_2 - a_2 b_1 + a_3 b_4 -a_4 b_3)^2  \\\\\n",
    "  &+& (a_1 b_3 - a_2 b_4 - a_3 b_1 +a_4 b_2)^2  \\\\\n",
    "  &+& (a_1 b_4 + a_2 b_3 - a_3 b_2 -a_4 b_1)^2.\n",
    "\\end{eqnarray*}\n",
    "\n",
    "The following code checks this; it can also be found in the test suite for the package.\n",
    "\n",
    "We first define the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "make_variable (generic function with 1 method)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define the variables α₁, ..., α₄, β₁, ..., β₄\n",
    "make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "variable_names = [make_variable(\"α\", i) for i in 1:4];\n",
    "\n",
    "append!(variable_names, [make_variable(\"β\", i) for i in 1:4]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create the Taylor objects (order 4, numvars=8)\n",
    "a1, a2, a3, a4, b1, b2, b3, b4 = set_variables(variable_names, order=4);    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 α₁ + 𝒪(‖x‖⁵)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 β₁ + 𝒪(‖x‖⁵)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute each term appearing in the above equation, and compare them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 α₁² β₁² + 1.0 α₁² β₂² + 1.0 α₁² β₃² + 1.0 α₁² β₄² + 1.0 α₂² β₁² + 1.0 α₂² β₂² + 1.0 α₂² β₃² + 1.0 α₂² β₄² + 1.0 α₃² β₁² + 1.0 α₃² β₂² + 1.0 α₃² β₃² + 1.0 α₃² β₄² + 1.0 α₄² β₁² + 1.0 α₄² β₂² + 1.0 α₄² β₃² + 1.0 α₄² β₄² + 𝒪(‖x‖⁵)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# left-hand side\n",
    "lhs1 = a1^2 + a2^2 + a3^2 + a4^2;\n",
    "lhs2 = b1^2 + b2^2 + b3^2 + b4^2;\n",
    "lhs = lhs1 * lhs2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 α₁² β₁² + 1.0 α₁² β₂² + 1.0 α₁² β₃² + 1.0 α₁² β₄² + 1.0 α₂² β₁² + 1.0 α₂² β₂² + 1.0 α₂² β₃² + 1.0 α₂² β₄² + 1.0 α₃² β₁² + 1.0 α₃² β₂² + 1.0 α₃² β₃² + 1.0 α₃² β₄² + 1.0 α₄² β₁² + 1.0 α₄² β₂² + 1.0 α₄² β₃² + 1.0 α₄² β₄² + 𝒪(‖x‖⁵)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# right-hand side\n",
    "rhs1 = (a1*b1 + a2*b2 + a3*b3 + a4*b4)^2;\n",
    "rhs2 = (a1*b2 - a2*b1 + a3*b4 - a4*b3)^2;\n",
    "rhs3 = (a1*b3 - a2*b4 - a3*b1 + a4*b2)^2;\n",
    "rhs4 = (a1*b4 + a2*b3 - a3*b2 - a4*b1)^2; \n",
    "rhs = rhs1 + rhs2 + rhs3 + rhs4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lhs == rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The identity is thus satisfied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Fateman test\n",
    "\n",
    "Richard J. Fateman, from Berkley, proposed as a stringent test\n",
    "of polynomial multiplication\n",
    "the evaluation of $s \\cdot (s+1)$, where $s := (1+x+y+z+w)^{20}$. This is\n",
    "implemented in\n",
    "the function `fateman1` below. We also evaluate the alternative form $s^2+s$ in `fateman2`,\n",
    "which involves fewer operations (and makes a fairer comparison to what\n",
    "Mathematica does). Below we have used Julia v0.4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "set_variables(\"x\", numvars=4, order=40);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fateman1 (generic function with 1 method)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fateman1(degree::Int)\n",
    "    T = Int128\n",
    "    oneH = HomogeneousPolynomial(one(T), 0)\n",
    "    # s = 1 + x + y + z + w\n",
    "    s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )\n",
    "    s = s^degree  \n",
    "    # s is converted to order 2*ndeg\n",
    "    s = TaylorN(s, 2*degree)\n",
    "    s * ( s+TaylorN(oneH, 2*degree) )\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.205706 seconds (198.63 k allocations: 8.903 MB, 2.13% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time f1 = fateman1(0) #for compilation; "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.694599 seconds (4.12 k allocations: 25.513 MB, 0.16% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time f1 = fateman1(20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another implementation of the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fateman2 (generic function with 1 method)"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fateman2(degree::Int)\n",
    "    T = Int128\n",
    "    oneH = HomogeneousPolynomial(one(T), 0)\n",
    "    # s = 1 + x + y + z + w\n",
    "    s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )\n",
    "    s = s^degree\n",
    "    # s is converted to order 2*ndeg\n",
    "    s = TaylorN(s, 2*degree)\n",
    "    return s^2 + s\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.006899 seconds (7.60 k allocations: 333.759 KB)\n"
     ]
    }
   ],
   "source": [
    "@time f2 = fateman2(0) #for compilation;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.291331 seconds (4.05 k allocations: 24.573 MB, 0.24% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time f2 = fateman2(20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "128358585324486316800"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getcoeff(f2,[1,6,7,20]) # coef x^1 y^6 z^7 w^{20}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "attributes": {
     "classes": [
      "julia"
     ],
     "id": ""
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "135751"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(TaylorSeries.size_table) # number of distinct monomials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tests above show the necessity of using `Int128`, that\n",
    "`fateman2` is about twice as fast as `fateman1`, and that the series has 135751 monomials on 4 variables.\n",
    "\n",
    "Mathematica version 10.3.0 requires 7.606437 seconds. Then,\n",
    "with TaylorSeries v0.3.0, our implementation of `fateman1` is about a factor ~2.7 faster,\n",
    "and `fateman2` is a factor 5.8 faster. (The original test by Fateman\n",
    "corresponds to `fateman1`, which avoids optimizations in `^2`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.8-pre",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.8"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
